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Abstract. We present a general model allowing one to 
calculate the distribution function of energetic particles 
(EPs) in the interstellar medium (ISM), and hence any rel- 
evant nuclear reaction rate, for any given time-dependent 
injection function, as well as in situations when the prop- 
agation coefficients themselves are non stationary. We re- 
view and provide physical interpretation of general formal 
solutions of a propagation equation taking into account 
energy losses, nuclear destruction and escape. Both one- 
zone and extended models are investigated. Our main goal 
is to provide standard and general tools for subsequent use 
in nuclear astrophysics, notably to calculate the gamma- 
ray emission and spallation rates associated with ener- 
getic events and specific acceleration mechanisms. Con- 
sidering the stationary limit of the model, we show that 
only time-dependent calculations can probe the density 
dependence of the physical processes under study, while 
steady-state models can in principle only give informa- 
tion about density-independent processes. 
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1. Introduction 

The energetic particles (EPs) present in the Galaxy can be 
probed through their electromagnetic and nuclear interac- 
tions within the interstellar medium (ISM). These inter- 
actions can lead to different types of radiation, lines (e.g. 
from nuclear de-excitation or pair annihilation) or contin- 
uum (e.g. synchrotron, brcmsstrahlung, inverse Compton 
or 7T° disintegration), as well as to the synthesis of sec- 
ondary nuclei such as Li, Be and B (LiBeB), whose as- 
trophysical significance is well established (e.g. Reeves et 
al. 1970; Meneguzzi et al. 1971; Fields et al. 1994; Reeves 
1994; Ramaty et al. 1996, Vangioni-Flam et al. 1996). 

Basically, the calculation of spallation and gamma-ray 
production rates at any place in the Galaxy requires the 
knowledge of i) the cross sections for relevant physical pro- 
cesses and ii) the distribution function of all the interact- 
ing EP species, i.e. their energy spectrum and number 



density. In this paper, we indicate how to calculate these 
distribution functions in some general situations where an 
exact formal solution exists. The different reaction rates 
can then be obtained by merely integrating the corre- 
sponding cross sections over the EP spectral densities, 
given the local ISM (target) density and chemical com- 
position. The knowledge of the local magnetic field and 
radiation background further allows one to compute the 
synchrotron and inverse Compton continuum emission in 
much the same way. 

In most previous nuclear astrophysics calculations, 
a steady-state assumption has been used, in which the 
steady injection of EPs in the region under study (gen- 
erally as the result of some particle acceleration process) 
is counterbalanced by a drift in energy space due mainly 
to Coulombian energy losses. The resulting EP distribu- 
tion functions are thus independent of time, as are conse- 
quently all the reaction rates. In this paper, we investigate 
the case when i) the injection of EPs and ii) the conditions 
of their propagation in the surrounding medium are func- 
tions of time. The former type of time dependence arises 
naturally when the acceleration conditions in a specific re- 
gion are changing, for instance when a supernova explodes, 
allowing for acceleration locally in space and time. As for 
the latter, it can happen when the density and/or com- 
position of the ambient medium varies, as for example in 
an expanding medium such as the interior of a supernova 
remnant (Parizot & Drury 1999a, b). 

In general, a steady-state model is sufficient provided 
that the time scale of the changes in the injection and/or 
propagation conditions is larger than the time needed for 
the EP distribution function to reach its equilibrium, that 
is, basically, the time scale of energy losses or escape out 
of the region considered. In this case, the history of EP- 
ISM interactions is fairly well reproduced by juxtaposing 
different phases in which the steady-state approximation 
holds. Such a strategy has been used extensively for in- 
stance to describe the Galactic evolution of light element 
abundances, as induced by GCR interactions with the ISM 
(most recently, Fields & Olive 1999; Vangioni-Flam et al. 
1999). Indeed, the spallation rates at time t depend on the 
current ambient metallicity and the current flux of GCR, 
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which are both functions of time, but with time scales pre- 
sumably larger than the GCR confinement time. There are 
some situations, however, where the steady-state approx- 
imation does not hold (Parizot et al. 1997a, b,c; Parizot & 
Drury 1999a, b), and a time dependent model is needed. 

It has to be noted that formal solutions of the propa- 
gation equation have been known for several decades (e.g. 
Syrovatskii, 1959; Jones, 1970; Reames, 1974; Ramaty, 
1974) and used extensively in the cosmic-ray transport 
theory, with additional features and complications, such 
as convection, re-acceleration, spatially dependent diffu- 
sion, etc. (Lerche and Schlickeiser 1981, 1988; Schlickeiser 
1986; Wang and Schlickeiser, 1987). However, the case 
when the energy losses and/or the diffusion coefficient are 
cxplicitely depending on time, and not only on space, has 
not yet been considered, to the best of our knowledge, 
while it is relevant to many applications in nuclear astro- 
physics (e.g. Parizot & Drury 1999a, b). Moreover, most of 
the situations of interest for spallative nucleosynthesis and 
gamma-ray line astronomy do not require the refinements 
mentioned, so we shall try to summarize and present in 
a clear and unified way the only results which one needs 
to develop a general time-dependent code for nuclear as- 
trophysics calculations. In particular, we shall emphasize 
the simple and intuitive physical meaning of the solutions, 
describing with special care the more original features, 
especially useful for numerical implementation, such as 
the temporal connection of solutions obtained in different 
phases of the EP propagation. 

2. Basic physical ingredients 

The general structure of a spallation or gamma-ray pro- 
duction calculation can be divided into three stages : ac- 
celeration (production of the EPs), propagation (trans- 
port in the surrounding medium) and interaction (with 
the surrounding matter or radiation field). These three 
logical stages are not necessarily separated in time, as (re- 
) acceleration may occur while the EPs are propagating, 
and interactions with the ambient medium just cannot 
be avoided at any stage of the process. In many situa- 
tions, however, the acceleration time scale of the particles 
is much shorter than their interaction and energy loss time 
scales, so that acceleration can actually be treated simply 
as injection (of EPs). This means that the acceleration 
itself is not affected by the conditions of propagation, and 
therefore the two stages are disconnected. In this case, we 
are solely concerned with the future of the EPs once they 
have been accelerated (or advected from a distant accel- 
eration site) and injected into the region of space under 
study. 

Accordingly, the first physical ingredient that we need 
to know is the injection function, qi(E,r,t), defined as 
the number of particles of species i introduced (injected) 
at energy E, location r and time t, per unit energy, volume 
and time (in (MeV/n)~ 1 cm _3 s _1 ). This injection function 



can be either postulated, in order to reproduce specific 
astronomical observations, or calculated as the output of 
some known acceleration mechanism. 

The treatment of the two subsequent stages, propa- 
gation and interaction, is also made easier by artificially 
separating them in time, which can be legitimated in most 
cases as discussed below. In any case, the influence of the 
propagation of the EPs on their distribution function can 
be derived from the knowledge of : i) spatial diffusion prop- 
erties, ii) drift properties in the energy space, iii) the local 
production rate (as secondary nuclei) and iv) the 'catas- 
trophic' loss rate. By catastrophic loss we mean the out- 
right disappearance of the particle, either by escape out 
of the region under study, destruction in a nuclear reac- 
tion (in-flight spallation) , or radio-active decay in the case 
of unstable particles. This is mathematically described by 
the total loss time, r* ot (i?, r, t), which depends on the nu- 
clear species considered, i, and is a priori a function of 
energy, location and time. It includes the escape, nuclear 
destruction and decay times according to : 
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where j(E) is the Lorentz factor of the particle. 

Strictly speaking, the escape time, rf sc , is relevant only 
to one-zone models, when the r variable is meaningless. 
Indeed, in spatially extended models, the escape of parti- 
cles out of the region under study is taken automatically 
into account through diffusion, and should not be included 
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The nuclear destruction time, r 2 D , is obtained from the 
total inelastic cross sections &ij for a projectile i in a 
target of species j, as : 
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where v(E) is the velocity of the particle and nj(r,t) is 
the local number density of the target nuclei of species j, 
at time t. 

Similarly, the production rate of the EP species k as 
secondary nuclei, through all the nuclear reactions i + j — > 
k with cross section aij-k{E), is given by : 

q' k (E,r,t) = y2 / 6E , n i {E',r,t)n j {r,t)<T i& , k {E')vi{E'), 
o Jo 

(3) 



where rn(E,r,t) is the distribution function of energetic 
nuclei of species i, i.e. their number per unit volume and 
energy (in (MeV/n) _1 cm -3 ). Although q^E, r,t) may be 
thought of as part of the injection function qi(E,r,t), it 
should be noted that it is actually the part of the injection 
rate which depends on the EP content of the ISM. This 
makes q 1 , mathematically more difficult to implement than 
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the primary injection function, qi, which is the same what- 
ever the local and current distribution functions may 
be. 

Apart from being created (synthesized) and removed 
by catastrophic losses, EPs are subject to energy losses 
which modify the shape of their distribution function. 
Whether coulombian, synchrotron, inverse Compton, adi- 
abatic or of any kind, these energy losses can gener- 
ally be expressed by means of an energy loss function, 
dE/dt = Ei[E, r, t), giving the energy loss rate (E l < 0) 
of nuclei i (in (MeV/n)s _1 ) as a function of energy, lo- 
cation and time. Note that the energy loss function can 
also include positive terms corresponding to (non diffu- 
sive) steady acceleration. 

The last physical ingredient that we need to know 
is the diffusion coefficient Di(E,r,t). Together with 
qi(E,r,t), r*; ot (E,r,t) and Ei(E,r,t), it allows one to 
calculate the distribution function of each EP species, 
Ui(E, r, t), at least in principle. The rate of any given reac- 
tion at any place and time is then obtained by integrating 
the relevant cross section over Ui(E,r,t) in the energy 
space, just as in Eq. (|3|), where the index k may represent 
any secondary nucleus or photon species (e.g. from a given 
gamma-ray line) in which we are interested. 

The above ingredients being assumed known, we can 
now write down the propagation equation and its solution 
in a few general cases. 

3. One-zone models with stationary propagation 
conditions 

In a number of cases, it is useless to study the spatial dif- 
fusion of the EPs, not because it is negligible, but because 
we are interested in the gamma-ray emission or secondary 
nuclei production integrated over a sufficiently large vol- 
ume to be thought of as a box without internal structure. 
A 'one-zone model' can then be used, in which the space 
coordinates are purely and simply dropped, and the EP 
distribution function, Ui(E, r, t), is replaced by its integral 
over the volume, namely the spectral density, Ni(E, t), de- 
fined as the differential number of EPs at energy E and 
time t (in (MeV/n) -1 ) in the whole box. For example, 
such models are appropriate for the evaluation of reaction 
rates in a homogeneous medium, or more generally in an 
instrument ally unresolved region of space. 

In a one-zone model, we have to solve the EP propa- 
gation equation in the following form : 



-^N i {E,t) + —{E i {E)N i {E,t)) 
= Q i (E,t) + Q' i (E,t)- 



Ni{E,t) 



(4) 



which merely expresses that the rate of change of the spec- 
tral density is equal to what is injected minus what is lost. 
Qi(E,t) and Q'^E^t) are the integrals of qi(E,r,t) and 
q[{E, r,t) (i.e. Eq. (g)) over the volume of the box, and 



the second term in the left hand side describes the drift 
in the energy space, due to energy losses. 

Let us now distinguish between primary EPs, for 
which Qi(E,t) 3> Q^Ejt), and secondary EPs, for which 
Qi(E,t) <C Qi(E,t). Primary EPs are typically protons, 
alpha particles and abundant metals such as C, N, O, 
Mg, Si, Fe, etc., whose production rates by spallation are 
usually much smaller than their injection rate from the 
acceleration of ambient materials. Energetic electrons are 
also mainly primary EPs. On the contrary, nuclei such 
as Li, Be and B, as most of the usual spallation prod- 
ucts, are secondary EPs. This is also the case of positrons 
and anti-protons, both of which may be important trac- 
ers of cosmic-ray propagation. Such a distinction simplifies 
the calculation both conceptually and mathematically. In- 
deed, we shall first solve the propagation equation for the 
primary nuclei, with the term Q'^E^t) dropped, and then 
use the spectral densities obtained in this ways to calculate 
the production rates of secondary nuclei. For EP species 
with Qi(E,t) and Q^(i?,i) of the same order of magni- 
tude, one has to solve first the propagation equation for 
the parent nuclei, deduce Q'^E^t) and then solve Eq. ^ 
again, in a (hopefully quickly convergent) iterative way. 

The formal solution of Eq. (0) reads : 
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where the function Ti(Eo,E) has an important physical 
meaning. It is defined for each EP species i as : 



Ti(Eo,E) 



dE' 



(G) 



'E Me'Y 

so it has the dimensions of a time, and describes the drift 
in energy of the EPs, due to the various energy loss mech- 
anisms, operating at a rate En. We shall call r the energy 
drift function, or drift time. It is clear from Eq. (|^) that 
Tj(i?o, E) is nothing but the time needed by a particle of 
species i to slow down from energy E to energy E. This 
makes the formal solution, Eq. (|J), very easy to under- 
stand. For a particle having a given energy, E, at time t, 
the drift function sets a one-to-one relation between the 
time at which it has been injected in the box, and the 
energy it then had. Eq. (^) thus merely says that the par- 
ticles found at time t at energy E are the collection (the 
mathematical sum) of all the particles injected at a higher 
energy, Eq > E, but at an earlier time, to, such that since 
this time, their energy has dropped exactly to energy E, 
as the result of the various energy loss mechanisms. The 
argument to = t — ti{Eo,E) can thus be interpreted as 
a mere retarded time, expressing the delay between the 
injection at (Eo,to) and the collection at (E,t) . 

As for the exponential factor in Eq. (g), it is equal to 
one in the absence of catastrophic losses (r to t - > oo), and 



4 



Etienne Parizot: Time dependent models for EP interactions in the ISM 



otherwise corrects the contribution of the EPs injected at 
energies Eq > E by weighting it by the probability of sur- 
vival during the time needed to slow down from energy 
Eq to energy E, that is during the drift time ti{Eo,E). 
This survival probability is Vi — exp(— (rj/Vtot))) where 
the average accounts for the energy dependence of catas- 
trophic losses, Ttot(E'), which are indeed varying during 
the energy drift. The time spent by the particles between 
energies E' and E' + dE' being dt(E') = dE'/E z (E'), one 
gets : 



Ttot 



dE' 



E Ei(E')r totii (E')- 



(7) 



which justifies the form of the exponential factor in 
Eq. (§). 

4. One-zone models with time-dependent 
conditions of propagation 

In the previous section we assumed that the energy loss 
rate, Ei(E), depended on the EP energy, but not on time. 
Now let us consider the case of ionization losses in a 
medium of density po, whose chemical composition is given 
by the relative abundance (by number) of each nuclear 
species, j: otj. Let dE/dx\j be the particle energy loss per 
unit grammage passed through in a medium of pure j nu- 
clei (expressed in g/cm 2 ). Then the energy loss rate of 
particles i in the ISM is given by : 



Mm = £ 

dx 



x aijpov. 



(8) 



As a consequence, the EP energy loss rate cannot be con- 
sidered constant whenever the density or composition of 
the propagation medium are functions of time. Now such 
situations are not unusual in astrophysics. One only has 
to think of the medium surrounding an active massive 
star: it is composed of successive layers of different den- 
sity and composition, resulting from successive phases of 
stellar wind (main sequence, red supergiant, Wolf-Rayet 
of type N, C, etc.). Such a situation has been considered 
by Parizot et al. (1997a, b,c). Likewise, EPs interacting in 
the interior of a supernova remnant will experience varia- 
tions in the ambient density and composition (due to the 
expansion of the remnant and the dilution of the ejecta 
by the swept-up material). This situation is addressed in 
detail in Parizot & Drury (1999a, b). 

To deal with this non-stationarity of the conditions of 
propagation, we divide the process into successive phases 
during which the density and composition of the tar- 
get can legitimately be considered constant. The solu- 
tion of the propagation equation restricted to any in- 
dividual phase is obtained in exactly the same way as 
above (Eq. (||)), and the only technical difficulty lies in 
the proper connection between the successive stationary 
solutions. To see how this works, let us consider the situ- 
ation of Fig. [|. 



We are interested in the EP distribution function at a 
succession of times t, during three distinct phases follow- 
ing the onset of the EP injection process, assumed to occur 
at t — (by definition, the ambient density and composi- 
tion are constant during each of these phases). The set of 
these distribution functions would typically allow one to 
calculate observables such as the gamma-ray line emission 
rate as a function of time. Let us first calculate the distri- 
bution function during phase 1, i.e for < t < t\, where 
t\ marks the end of the first phase. Since the conditions of 
propagation are constant during each phase, the EP dis- 
tribution function at any time before t\ is readily given 
by Eq. (|^). Dropping the species index, i, for convenience 
and noting Ei(E) and ti(Eo,E) the energy loss function 
and the energy drift function during phase 1, we get : 



N(E,t) 



1 



\Ei(E)\ 



/ dE'Q(E',t-T 1 (E',B)yp 1 {E',E), 

J E 



(9) 



where Vi(E',E) = cxp(— (ri/rtot,i)) is the probability 
(calculated under the conditions of phase 1) for an EP 
injected at energy E' to reach energy E without escaping 
or being altered by nuclear reactions (see Eq. (Q)). 

Note that since the injection of EPs began at time 
t = 0, the injection rate is null at negative times, so that 
Q(E',t — Ti(E' , E)), considered as a function of the in- 
jection energy, E 1 , vanishes for all energies greater than a 
maximum value, Eq, defined by : 



n(E ,E)=t. 



(10) 



Indeed, any particle injected at a higher energy would need 
a time greater than t to slow down to energy E, and con- 
sequently cannot contribute to N(E,t). The upper limit 
in the integral of Eq. (g) can thus be replaced by the en- 
ergy E (E,t), implicitly defined in Eq. (|l0|). In order to 
shorten the writing of equations, we define the retarded 
injection function, fi(E',E), as : 

f 1 (E',E)=Q(E',t-T 1 (E',E))xV 1 (E l ,E), (11) 

whose physical meaning is straightforward, so that Eq. <JqJ> 
can be re-written as : 



N(E,t) = 



1 



E^E)] 



E (E,t) 



h{E',E) dE'. 



(12) 



Let us now calculate the EP spectral density, N(E,t), 
during phase 2, i.e. at times t\ < t < t 2 - One can still 
define a limiting energy, Ei(E,t), such that the particles 
having this energy at time t\ (end of phase 1, beginning 
of phase 2) have exactly energy E at time t. Ei(E,t) is 
defined similarly to Eq in Eq. (|l^) , except that the energy 
drift function, r, now corresponds to the specific propaga- 
tion conditions prevailing during phase 2 : 



T 2 (E 1 ,E)=t-t 1 . 



(13) 
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Fig. 1. Diagram corresponding to an EP propagation process divided into three phases, during which the ambient 
density and chemical composition, and thus the energy loss rate, are assumed constant. Instants t\ and t 2 correspond 
to the changes of phase, and t is the 'current' time, at which we want to calculate the EP energy spectrum. The energy 
axis goes backward in time, because the energy losses produce a continual decrease of the energy of particles. Particles 
having energy E at time t (whose density, N(E,t), we are looking for), had energy E 2 at time t 2 , E\ at time ti, and 
Eo at time to — 0. These limiting energies of course depend on E and t. The diagram can be easily generalized to 
as much phases as one wishes. The contribution of all phases to the current spectral density must be added, that of 
phase i being limited to a domain of energy stretching from to Ei (see text). 



As above, any particle injected at an energy E' > E\ 
cannot contribute to the spectral density N(E, t) unless 
it has been injected earlier than t\, i.e. during phase 1. 

Had the injection of EPs began at time ti, we could 
easily write the solution of the propagation equation dur- 
ing phase 2, just as above : 



N(E,t) 



1 



\E 2 (E)\ J E 
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It is given by Eq. (or its condensed version Eq. ( |l2| 

1 



N{E X {E, t)M 



\Ei{E x (E,t))\ 

fE {E,t) y lb > 

x / f 1 (E',E 1 (E,t))dE', 

J Ei_{E,t) 



where we recognized that 

E (Ei(E,t),t 1 )=E Q {E,t). 



(17) 



where the retarded injection function is now 

f 2 (E', E) = Q(E>, t - t 2 (E', E)) x V 2 (E', E), (15) 

and t 2 (E',E) and V 2 (E' \E) are the new energy drift 
function and EP survival probability, acknowledging the 
new target density and composition (i.e. corresponding to 
phase 2). But to obtain the actual EP spectral density, we 
still have to add the contribution of all the EPs injected 
during phase 1. Now by definition of Ei(E, t), N(E, t) col- 
lects at time t all the EPs which had energy E\ (E, t) at 
the end of phase 1, i.e. at time t\. But their number per 
unit energy, N(Ei(E,t),ti), has already been calculated. 



Before adding this contribution from phase 1 to the 
spectral density of the EPs at time t, we still need 
to take into account the catastrophic losses of particles 
since the end of phase 1. This is achieved by multiply- 
ing N(Ei(E,t),t\) by the survival probability (under the 
conditions of phase 2) from energy Ei(E,t) to energy E. 
Finally, we must remember that the distribution function 
N(E,t) gives the differential number of EPs at energy 
E, that is the number of EPs between E and E + dE. 
Accordingly, we must re-scale N(Ei(E),ti) by a factor 
\E2(Ei(E, t))\/\E2(E)\ representing the contraction (or 
dilation) of the 'comoving' energy interval. 

Summing up all this 'book keeping', we finally write 
the spectral density of the EPs at any time t during 
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phase 2 as : 



N(E,t) 



\E 2 (E)\ 



E 1 {E,t) 



f 2 (E',E)dE> 



N(E 1 {E,t),t 1 )P 2 {E 1 (E,t),E) 



\MEi(E,t))\ 



\E 2 (E)\ 



(18) 



or, substituting from Eq. (|T^) and conveniently abbrevi- 
ating Ei(E, t) to Ei, 



N(E,t) = -j 



\E 2 (E)\ 
\E 2 (Ei 



[ 1 f 2 (E',E)dE> 
Je 



(19) 



r ,^2(4£) Ei)dE' 
\Ei(Ei)\ J Ex 



Note that we could have alternately obtained Eq. ( |18[ ) 
by integrating the propagation equation, (Q), beginning at 
time t\ and replacing the injection function Q(E, t) by an 
effective injection function, Q c ff, defined by : 



Q cS (E, t) = Q(E, t) + N(E, h)6{t - h). 



(20) 



Although mathematically equivalent, we feel however that 
the derivation given above gives a deeper physical insight 
of the situation. 

We finish the study of the multiphase case by comput- 
ing the EP spectral density during phase 3. The calcula- 
tion leading to Eq. (|l^) generalizes straightforwardly if we 
introduce the energy, E 2 (E, t) , that a particle of current 
energy E (i.e. at time t) had at the time t 2 marking the 
end of phase 2 and the beginning of phase 3. As above, 
E 2 (E,i) is defined by : 



t 3 (E 2 , E) = t - t 2 , 



(21) 



where the subscript '3' refers to the conditions of propa- 
gation during phase 3. 

Recursively adding the contributions of each phase to 
the spectral density at time t, we obtain : 

fE 2 

\E 3 (E)\N{E,t) = / h(E',E)&E' 

JE 

+ 7M£ 2 ,£) / 1 f 2 (E',E 2 )dE' 

\E 2 (E 2 ) | Je 2 



\E 2 (E 2 )\ 



< I h{E',E{)&E', (22) 

'-Ei 



with straightforward generalization to any number of 
phases. 



5. Extended models with homogeneous spatial 
diffusion 

Although one-zone models are both powerful and easy to 
handle, gamma-ray observations with high angular reso- 
lution might make it also necessary to use extended mod- 
els, allowing one to study inhomogeneous astrophysical 
sites. As an example, one may be interested in situations 
where energetic particles are accelerated in a small region 
of space and diffuse in the surrounding medium, exploring 
sites with various densities and magnetic properties. The 
EP distribution function, and hence the gamma-ray emis- 
sion, should thus be different in neighboring regions which 
can be resolved by the instruments. Moreover, as was dis- 
cussed in Parizot (1997), the detailed study of spatial dif- 
fusion is sometimes indispensable, even in homogeneous 
media, if the loss of EPs by escape out of the confinement 
region has to be taken into account. 

Let us then write the propagation equation of the EPs 
in the case of an extended model. We are now interested 
in the EP distribution function, m(E,r,t), rather than 
in the integrated spectral density, Ni(E,t). It satisfies an 
equation similar to Eq. (Q) , with new terms describing the 
spatial transport of particles. In a number of applications, 
convection can be neglected, because one can work in the 
referential of the interstellar plasma, or average over large 
enough regions of space for the plasma to be considered 
globally at rest. Keeping only the diffusion term, we then 
have the EP propagation equation in the extended case : 

^n,(r, E, t) + ^(Eifflmir, E, t)) = qi (r, E, t) 



m(r,E,t) 



Ttot 



V(D(r,E,t)Wm{r,E,t)). (23) 



This equation can of course be solved numerically, but 
it is much better (be it only regarding computation time) 
to solve it formally once for all, for any injection function 
<7i(r, E, i), as in the one-zone case. To do this, however, we 
shall assume that the diffusion coefficient does actually not 
depend on the location, r, and time, t. On the other hand, 
it can still depend on (and be any function of) the energy, 
as has to be the case in most realistic astrophysical situa- 
tions. Note that although the condition D(E, r, t) = D(E) 
may be too restrictive in some specific occasions, there 
will be a way around it in most cases. Indeed the diffusion 
conditions (gathered mathematically in the diffusion co- 
efficient, D) do not generally change smoothly, but rather 
sharply, while passing from one type of environment to 
another, like from the ISM to a dense cloud or from a 
superbubble to the ISM. In such cases, one can use a pro- 
cedure analogous to that of Sect. [I| and solve the propaga- 
tion equation locally, in regions with homogeneous diffu- 
sion, and connect the solutions at the boundaries in much 
the same way as we did above in the multi-phase case. In 
practical terms, this amounts to adding to the injection 
function, qi(E,r,t), a term of the form J s 5(r — r )dS, 
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corresponding to the flux of particles across the boundary 
of each region. This flux 'inherited' from the propagation 
of EPs in neighboring regions is exactly analogous to the 
spectral density 'inherited' from the propagation during 
the earlier phase, as we encountered it in Sect. [|. 

The solution of Eq. (^) (with homogeneous diffusion) 
can be given in a convenient way by introducing a new 
function, Xi(Eo, E), in the spirit of the energy drift func- 
tion given by Eq. (0) : 



Je E(E' 



(24) 



Physically, Xi{E ,E) represents the effective value of l DV 
when the particle slows down from energy Eq (at which it 
has been injected) to energy E (which it has at time t). 
This is easily seen by noting that dE' / E(E') represents 
the time passed by the particle between energies E' + 
dE' and E'. The introduction of this 'effective' diffusion 
parameter is necessary because the energy of the particle 
decreases with time, while the diffusion coefficient depends 
on energy. In the case when D(E) = Dq is constant, we 
have in fact Xi(Eo,E) — D Ti(E , E). Note finally that 
Xi(Eo, E) has the dimensions of a surface, so we shall call 
it the diffusion section of the particle (when passing from 
E to E). 

We can now give the EP distribution function at en- 
ergy E, location r and time t, in the case of a homoge- 
neous, though energy dependent diffusion coefficient : 

1 



nAE, r.t) = — 

\Ek(E)\ 

dE' J dr' Mr';E',E)g t (r',r;E',E), (25) 

where fi is the 3D generalization of the retarded injection 
function (Eq. (O)) : 



Mr'; E', E) = qi (E', r' , t - n {E', E)) 



x exp 



dE" 



E' Ei(E")Tt ot ,i(E") 



). (26) 



and §i is given by : 
gi (r',r;E',E) 



8(7T X (E' 7 E)r/' 



■ exp 



\\ r -r II 
i X (E',E) 



(27) 



This last function proves to have a very simple physi- 
cal interpretation. To see this, let us consider the simplest 
equation of diffusion, for identical particles in a homoge- 
neous and isotropic medium, dN/dt — DAN — 0, where 
N is the number density of particles, D is the diffusion 
coefficient, and A denotes the Laplacian operator. In an 
infinite medium, if No particles are initially injected at 
r = ri n , then we have the well-known solution : 

II™ m 112- 



N(r,t) 



N 



8(7rDi)3/2 



exp 



ADt 



(28) 



One way to look at this solution, emphasizing the stochas- 
ticity of the diffusion process, is to say that if a particle 
is introduced in r ln at t = 0, the density of probability 
that it be in r at time t is Vt{r m ,r) — N{r,t)/No- Now 
comparing the expressions for Vt(ri n , r) and g(r', r; E' , E) 
(given by Eq. (p7j)), while remembering that x(E',E) is 
the effective 'Dt' during the time it takes the particles to 
slow down from E' to E, it is clear that 



gi (r',r;E',E) = V TiE , tE) (r',r) 



(29) 



or, to say it in words, the function <?i(r', r; E', E) appear- 
ing in Eq. ( p5| ) is nothing but the probability for a particle 
injected at location r' at energy E' to be found in r at 
energy E. 

With the definition hi — fi x g iy we can now re- write 
the solution ( |25| ) of the propagation equation in the sim- 
plest way : 

m(E,r,t) = — / dE' [ dr' hi(r',r:E',E), 
lK ' \Ei(E)\ J E J 

(30) 

where hi should be thought of as the 'full' retarded injec- 
tion function, i.e. retarded in both space and time (factors 
<?i and fi, respectively). The physical interpretation of the 
solution is then obvious: the particles which one finds at 
energy E, time t and location r, are the collection of all the 
particles which have been injected at a higher energy E' , 
anywhere in the volume considered (integration over r'), 
but at a time in the past such that, since their injection, 
the energy losses have brought their energy down from E' 
to exactly E. In addition, each of these contributions must 
be weighted by the probability that the particle survived 
during this time (exponential factor in Eq. (p6|)) and that, 
during this time again, it diffused from its injection site, 
r', to the place now probed, in r (function gi, Eq. (p7|)). 

6. The stationary limit of the models 

As was discussed in the introduction, some astrophysi- 
cal situations need not be studied in a time-dependent 
scheme. In steady-state, after a transitory regime lasting 
about the time-scale of energy losses, the spectral density 
of the EP species i reaches its equilibrium value, Ni(E), 
which satisfies : 

^-(E^NiiE)) = Qi(E) - (31) 

Otj Ttot 

Integrating this equation directly or, equivalently, taking 
the stationary limit of the general, time-dependent solu- 
tion given above (Eq. (||)), one obtains : 



Ni(E) 



f + Qi(E )V i (Eo,E)dE , (32) 
l-fM-frjl JE 



where the survival probability Vi has the same expression 
as above. 
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It is important to note that, in the steady-state approx- 
imation, the spectral density of stable energetic particles 
is inversely proportional to the ambient density, provided 
that the escape can be neglected, i.e. for low energy par- 
ticles or in the case of a thick target. Indeed, the Coulom- 
bian energy loss rate, Ei(E), is always proportional to the 
ambient density, no, while the catastrophic loss time, Ttot, 
is then equal to the nuclear destruction time, td, which 
according to Eq. (||), behaves as njj~ . As a result, the sur- 
vival probability, Vi, is independent of no, and inspection 
of Eq. ([32j) shows that Ni(E) scales as n^ 1 . Now the nu- 
clear reaction rates, due to the interaction of the EPs with 
the ISM, are obtained by integrating the cross-sections 
over the spectral density : 

r+oc 

Qk = J2 n j N i {E)a ij ^ k {E!)v i (E)dE, (33) 
i,i Jo 

where the notations are the same as in Eq. (||). 

It is clear, then, that the ii^ 1 dependence of Ni is 
exactly compensated by the factor rij, which is just no 
times the relative abundance of nuclei j in the ISM. As 
a conclusion, all the nuclear reaction rates are indepen- 
dent of density in the steady-state approximation. This is 
worth emphasizing, because it is one fundamental phe- 
nomcnological distinction between steady-state and time- 
dependent models. In essence, steady-state models cannot 
acknowledge density dependences. This is no longer true 
for time-dependent models, as we have shown elsewhere 
on various examples (Parizot et al. 1997a, b,c; Parizot & 
Drury 1999a,b). 

7. Conclusion 

In conclusion, we have reviewed the solutions of the dif- 
ferential equation describing the propagation of energetic 
particles (EPs) in the ISM, in the case of both one-zone 
and extended models, taking into account homogeneous 
spatial diffusion. The solutions obtained admit a sim- 
ple physical interpretation, which we emphasized intro- 
ducing intermediate functions : the energy drift function, 
Ti(E ,E), and the diffusion section, Xi(-E , E), which both 
depend on the nuclear species of the EPs. The first one 
describes the time it takes a particle to slow down from 
its injection energy, Eq, to its current energy, E, because 
of energy dependent energy losses. As for the latter, it 
plays the role of an effective 'diffusion coefficient-times- 
duration', taking into account the energy dependence of 
the diffusion coefficient. We have used these functions to 
construct retarded injection functions which sum up the 
energy history of the EPs and allow one to write the so- 
lution of the propagation equation in a simple and intelli- 
gible way. 

We also studied the case when the conditions of prop- 
agation vary, which arises when either the composition or 
the density of the ambient medium are not constant. We 



have shown how it is possible to use the standard results 
in this case, dividing the process into as much successive 
phases as required so that the conditions of propagation 
can be considered constant during each phase. We outlined 
the detailed procedure for fitting the phases altogether, re- 
lying on, and emphasizing again, the physical content of 
the mathematics involved. The procedure can be used to 
solve any type of time-dependent astrophysical situation 
involving the interaction of energetic particles in the inter- 
stellar medium. Indeed, the formal solutions derived here 
can be easily integrated to obtain the distribution function 
or spectral density of the EPs, from which the nuclear or 
electromagnetic reaction rates are deduced directly from 
the knowledge of the cross sections involved. 

Finally, we considered the stationary limit of the mod- 
els, stressing that the assumptions inherent in steady-state 
make it impossible in principle to investigate the influ- 
ence of the ambient density on the processes under study. 
This is one crucial superiority of the time-dependent mod- 
els which, for this reason, cannot be avoided in a great 
number of astrophysical situations. As an example, Pari- 
zot et al. (1997a, b,c) have studied the time-dependent 
gamma-ray line emission induced by the winds of mas- 
sive stars, and Parizot & Drury (1999a, b) calculated the 
spallative nucleosynthesis within an expanding supernova 
remnant, taking the dynamics of the process into account. 
In each case, the influence of the ambient density has been 
investigated, and in each case it proved very significant. 
We believe that the improvement of the spatial resolu- 
tion of gamma-ray detectors as well as the accumulation 
of data on variable gamma-ray sources will make the use 
of time-dependent models more and more necessary. 

In the same spirit, as the Galactic chemical evolution 
models are getting more and more precise, the question 
should be raised whether steady-state models are still ap- 
propriate in all the situations. In particular, recent studies 
of the evolution of the spallation products Li, Be and B, 
have focused mainly on the very early stages of Galactic 
evolution, when the metallicity of the ISM was less than 
one thousandth of the solar metallicity. Depending on the 
models, this corresponds to an age of the order of the con- 
finement or the energy loss time scales. As a consequence, 
one should be careful in assuming steady-state, because 
the transitory regime may have not yet totally decreased. 
In addition, the injection of cosmic rays is usually assumed 
to follow closely the supernova activity in the Galaxy, but 
it should be kept in mind that the supernova explosion 
rate may have changed on a very small time scale in the 
early Galaxy, which would be one more reason to be care- 
ful with steady-state models, and to investigate possible 
significant time-dependent effects. 
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